Hello everyone,

I am trying to reproduce some results from the paper:

Peristaltic flow of a magnetohydrodynamic nanofluid through a bifurcated channel

I want to generate:

• Fig 1 (channel geometry)

• Fig 8 (pressure gradient vs axial coordinate ξ)

• Fig 11 (streamlines)

• Table 1 (resistance values)

I am using Maple 2018.

Below is my Maple code:
p1_ravi_sir.mw

Problem:

The plots do not match the paper

p1_ravi_sir.pdf

Why does my resistance value not change with M

 Could someone please help me:

• Correct the Maple code

• Solve equations (18–22) properly in Maple

• Generate the correct plots 

Thank you very much.

Hi !

When I use the "expand" command on a summation, I've noticed that sometimes the command works and

sometimes it doesn't. I looked through the previous commands to see which one might be causing this problem.

I found that if I use the "changevar" command  from the "student" package,  (deprecated package,I know )

the "expand" command works correctly.     Strange behavior, isn't it ?       I have Maple 2018 on Windows 10.

Regards  !

There is a small  example in the file.

expand-problem.mw

The inscribed square problem, also known as the Toeplitz conjecture, is an unsolved quastion in geometry: Does every plane simple closed curve (Jordan curve) contain all four vertices of some square? This is true if the curve is convex or piecewise smooth and in other special cases. The problem was proposed by Otto Toeplitz in 1911. For detailes see  https://en.wikipedia.org/wiki/Inscribed_square_problem

The Inscribed_Square procedure finds numerically one or more solutions for a curve defined by parametric equations of its boundary or by the equation F(x,y)=0. The required parameter of procedure  L  is the list of equations of the boundary links or the equation  F(x,y)=0 . Optional parameters:  N  and  R . By default  N='onesolution' (the procedure finds one solution), if  N  is any symbol (for example  N='s'), then more solutions.  R  is the range for the length of the side of the square (by defalt  R=0.1..100 ).

The second procedure  Pic  visualizes the results obtained.

The codes of the procedures:

restart;
Inscribed_Square:=proc(L::{list(list),`=`},N::symbol:='onesolution',R::range:=0.1..100)
local D, n, c, L1, L2, L3, f, L0, i, j, k, m, A, B, C, P, M, eq1, eq2, eq3, eq4, eq5, eq6, eq7, eq8, eq9, sol, Sol;
uses LinearAlgebra;
if L::list then
L0:=map(p->`if`(type(p,listlist),[[p[1,1]+t*(p[2]-p[1])[1],p[1,2]+t*(p[2]-p[1])[2]],t=0..1],p), L);
c:=0;
n:=nops(L);
for i from 1 to n do
for j from i to n do
for k from j to n do
for m from k to n do
A:=convert(subs(t=t1,L0[i,1]),Vector): 
B:=convert(subs(t=t2,L0[j,1]),Vector):
C:=convert(subs(t=t3,L0[k,1]),Vector): 
D:=convert(subs(t=t4,L0[m,1]),Vector):
M:=<0,-1;1,0>;
eq1:=eval(C[1])=eval((B+M.(B-A))[1]);
eq2:=eval(C[2])=eval((B+M.(B-A))[2]);
eq3:=eval(D[1])=eval((C+M.(C-B))[1]);
eq4:=eval(D[2])=eval((C+M.(C-B))[2]);
eq5:=eval(DotProduct(B-A,B-A, conjugate=false))=d^2;
sol:=fsolve([eq1,eq2,eq3,eq4,eq5],{t1=op([2,2,1],L0[i])..op([2,2,2],L0[i]),t2=op([2,2,1],L0[j])..op([2,2,2],L0[j]),t3=op([2,2,1],L0[k])..op([2,2,2],L0[k]),t4=op([2,2,1],L0[m])..op([2,2,2],L0[m]),d=R});
if type(sol,set(`=`)) then if N='onesolution' then return convert~(eval([A,B,C,D],sol),list) else c:=c+1; Sol[c]:=convert~(eval([A,B,C,D],sol),list) fi;
 fi; 
od: od: od: od:
Sol:=fnormal(convert(Sol,list),7);
print(Sol);
ListTools:-Categorize((X,Y)->`and`(seq(is(convert(X,set)[i]=convert(Y,set)[i]),i=1..4)) , Sol);
return map(t->t[1],[%]);
else
A,B,C,D:=<x1,y1>,<x2,y2>,<x3,y3>,<x4,y4>:
M:=<0,-1;1,0>:
eq1:=eval(C[1])=eval((B+M.(B-A))[1]):
eq2:=eval(C[2])=eval((B+M.(B-A))[2]):
eq3:=eval(D[1])=eval((C+M.(C-B))[1]):
eq4:=eval(D[2])=eval((C+M.(C-B))[2]):
eq5:=eval(LinearAlgebra:-DotProduct((B-A,B-A), conjugate=false))=d^2:
eq6:=eval(L,[x=x1,y=y1]):
eq7:=eval(L,[x=x2,y=y2]):
eq8:=eval(L,[x=x3,y=y3]):
eq9:=eval(L,[x=x4,y=y4]):
sol:=fsolve({eq1,eq2,eq3,eq4,eq5,eq6,eq7,eq8,eq9},{seq([x||i=-2..2,y||i=-2..2][],i=1..4),d=R});
eval([[x1,y1],[x2,y2],[x3,y3],[x4,y4]], sol):
fi;
end proc:

Pic:=proc(L,Sol,R::range:=-20..20)
local P1, P2, P3, T;
uses plots, plottools;
P1:=`if`(L::list,seq(`if`(type(s,listlist),line(s[],color=blue, thickness=2),plot([s[1][],s[2]],color=blue, thickness=2)),s=L), implicitplot(L, x=R,y=R, color=blue, thickness=2, gridrefine=3));
P2:=polygon(Sol,color=yellow,thickness=0);
P3:=curve([Sol[],Sol[1]],color=red,thickness=3):
T:=textplot([[Sol[1][],"A"],[Sol[2][],"B"],[Sol[3][],"C"],[Sol[4][],"D"]], font=[times,18], align=[left,above]);
display(P1,P2,P3,T, scaling=constrained, size=[800,500], axes=none);
end proc:

Examples of use:

The curve consists of a semicircle, a segment and a semi-ellipse (find 1 solution):

L:=[[[cos(t),sin(t)],t=0..Pi],[[t,0],t=-1..0],[[0.5+0.5*cos(t),0.8*sin(t)],t=Pi..2*Pi]]:
Sol:=Inscribed_Square(L);
Pic(L,Sol);

The procedure finds 6 solutions for a non-convex pentagon:

 L:=[[[0,0],[9,0]],[[9,0],[8,5]],[[8,5],[5,3]],[[5,3],[0,4]],[[0,4],[0,0]]]:
Sol:=Inscribed_Square(L,'s');
plots:-display(Matrix(3,2,[seq(Pic(L,Sol[i]),i=1..6)]),size=[300,200]);

For an implicitly defined curve, only one solution can be found:

L:=abs(x)+2*abs(y)-sin((2*x-y))-cos(x+y)^2=3:
Sol:=Inscribed_Square(L);
Pic(L,Sol);

               
See more examples in the attached file.

Inscribed_Square.mw

The system below obviously has the unique solution  t=3*Pi/2 , but Maple doesn't return any solution. I wonder if this bug persists in recent versions of Maple?

solve({cos(t)=0,sin(t)=-1,t>=0,t<=2*Pi}, t);

   #  NULL

how to get the solution for this ode i am getting error 

in the pdf for the same parameter values getting plots but here it is showing error

santosh_sir_work.pdf

santosh_sir_work.mw

How to get a 2D plot for the following command: 
error_in_2d_plot_cf_and_nu.mw

Good day to you all ..

I built a routine to solve a machine-job allocation problem. This involves binary decision variables and an objective function to minimize the maximum completion time (makespan).The model assumes that each job, n, has a sequence of operations that must be processed in a specific order on a specific machine, m.

Consider (for example) the case of 5 machines and 3 jobs. Every job, 1, 2, and 3, visits machines 1 to 5 in order, with the processing times specified by the matrix, P. Maple successfully returns the required solutions for various combinations up to this scale.

However, when I encounter the 8-machine, 5-job problem, Maple stops (after 30 minutes of processing time) and returns the following error:

"Maximal depth, 182, of branch and bound search is too small; use depthlimit option to increase depth."

I have attempted to modify the depth limit (to a value exceeding 182) but without success ... I cannot resolve this problem. Can somebody advise on how to fix this - if possible? The worksheet in question is attached.

Thanks for reading.

MaplePrimes_Nov_16.mw

I am trying to plot a contour graph for my problem for (psi) function in the particular boundary, and even though it's working, but the contour  plot is not appearing at the end. Could anyone help me with the code to get proper graph in the specified boundary. 

i have ploted the graph in python i got a plot similar to that i am trying maple but i am not able to plot it. could any one help me to solve.

contour_plots_error_in_wavey_flow.mw

Good afternoon.

Here is (what appears to be) a simple assignment problem that involves scheduling 3 jobs through 2 machines, and I wish to minimize the total time required to complete all 3 jobs. Each job must be processed first on machine #1 before going to machine #2,and each machine can only handle one job at a time.

I have determined the optimal sequence to minimize the makespan (i.e., the completion time of the last job on machine #2) using a manual method and I was wondering if someone could review this (please see attached) and recommend a structured LP method or let me know if such a procedure exists in Maple's Optimization package.

Any help will be appreciated! 

Thanks for reading.

MaplePrimes_Nov_4.mw

Hi,
I want to have a colorbar aside the curve, but I can't. (My Maple version is 2018).

restart:
with(plots):

# --- constants ---
qc := 0.4:
qh := 0.64:
alpha := 0.8:

# --- parameter ranges ---
deltac_min := 0.01: deltac_max := 0.99:
Tch_min := 0.1:     Tch_max := 1.1:

# --- numerical safety ---
epsDen := 1e-12:
Mcap := 3.0: # cap to avoid infinity

# --- denominator ---
Den := (deltac, Tch) ->
       -2*Tch*deltac*qh - 2*Tch*deltac + 2*Tch*qh
       + 2*deltac*qc + 2*Tch + 2*deltac:

# --- safe functions ---
M1safe := (deltac, Tch) ->
    piecewise(Den(deltac,Tch) > epsDen,
              min(2/sqrt(Den(deltac,Tch)), Mcap),
              Mcap):

M0safe := (deltac, Tch) ->
    piecewise(Den(deltac,Tch) > epsDen,
              min(2*alpha/sqrt(Den(deltac,Tch)), Mcap),
              Mcap):

# --- density plots ---
p1 := densityplot(
       M1safe(deltac,Tch),
       deltac = deltac_min .. deltac_max,
       Tch   = Tch_min .. Tch_max,
       grid  = [200,200],
       style = patchnogrid,
       colorstyle = HUE,
       axes = boxed,
       labels = ["δ_c", "T_ch"],
       title = sprintf("M1(δ_c, T_ch)  (capped at %.2f)", Mcap)
     ):

p2 := densityplot(
       M0safe(deltac,Tch),
       deltac = deltac_min .. deltac_max,
       Tch   = Tch_min .. Tch_max,
       grid  = [200,200],
       style = patchnogrid,
       colorstyle = HUE,
       axes = boxed,
       labels = ["δ_c", "T_ch"],
       title = sprintf("M0(δ_c, T_ch)  (capped at %.2f)", Mcap)
     ):

# --- Contour Den = 0 (to mark physical boundary) ---
cont := contourplot(
          Den(deltac,Tch),
          deltac = deltac_min .. deltac_max,
          Tch   = Tch_min .. Tch_max,
          contours = [0],
          color = black,
          thickness = 2
        ):

# --- Overlay contour ---
p1_final := display(p1, cont):
p2_final := display(p2, cont):

# --- Display side by side ---
display(array([p1_final, p2_final]), scaling = constrained);

@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@

restart;
with(plots):

# === Constants ===
qc := 0.4:
qh := 0.64:
alpha := 0.8:

# === Denominator ===
Den := (deltac, Tch) ->
     -2*Tch*deltac*qh - 2*Tch*deltac + 2*Tch*qh
     + 2*deltac*qc + 2*Tch + 2*deltac:

# === M1 and M0 ===
M1 := (deltac, Tch) -> piecewise(Den(deltac,Tch)>0, 2/sqrt(Den(deltac,Tch)), undefined):
M0 := (deltac, Tch) -> piecewise(Den(deltac,Tch)>0, 2*alpha/sqrt(Den(deltac,Tch)), undefined):

# === Domain ===
deltac_min := 0.05:  deltac_max := 0.95:
Tch_min := 0.01:     Tch_max := 1.10:   # Physically realistic range

# === 3D Plots ===
p1 := plot3d(M1(deltac,Tch),
             deltac=deltac_min..deltac_max,
             Tch=Tch_min..Tch_max,
             grid=[80,80],
             axes=boxed,
             orientation=[-130,45],
             scaling=constrained,
             style=surfacecontour,
             color=green,
             labels=["δ_c", "T_ch", "M1"],
             labelfont=[TIMES,12],
             title="M1(δ_c, T_ch) for qc=0.4, qh=0.64"):

p2 := plot3d(M0(deltac,Tch),
             deltac=deltac_min..deltac_max,
             Tch=Tch_min..Tch_max,
             grid=[80,80],
             axes=boxed,
             orientation=[-130,45],
             scaling=constrained,
             style=surfacecontour,
             color=blue,
             labels=["δ_c", "T_ch", "M0"],
             labelfont=[TIMES,12],
             title="M0(δ_c, T_ch) for α=0.8, qc=0.4, qh=0.64"):

# === Display vertically (M1 above M0) ===
display(Array(1..2, [p1, p2]), insequence=true, scaling=constrained);
display(Array([[p1], [p2]]));

Good day.

Using the iterative map routine, the location of the bifurcation points of the logistic function can be determined using the plot (see attached).

I was wondering .. is it possible to estimate and output the locations of these points and the range of the function? I would like to explore the bifurcation behavior for the modified function so, this would be a great help. 

In this case, the location of the points are: (1.00,0.00), (3.00, 0.66), (3.45, 0.44), (3.45, 0.85) and the range is [0,4]. 

Thanks for reading!

MaplePrimes_Oct_16.mw

Hello everyone,

I am working on reproducing Fig. 12 (Isotherms) and Fig. 13 (Streamlines) from the attached paper, cacity_paper_work.pdf which deals with natural convection of GO–MgO/silicone oil hybrid nanofluid inside square and H-shaped cavities under a periodic magnetic field.

Equations (9–11) in the paper describe the steady-state 2D flow and heat transfer in terms of u(x,y), v(x,y), and θ(x,y). I have already set up the full PDE system, parameters, and boundary conditions (BCs) in Maple.

Now, iusing pdsolve(..., numeric, method=fd), 

Solve these equations for both square and H-shaped cavities using the given ICs and BCs, and

Generate Isotherm and Streamline plots (similar to Figs. 12 and 13 in the paper).

I have attached:

The PDF (showing equations and figures), and

My Maple worksheet containing the PDE definitions and parameters.

I would really appreciate guidance on:

Whether dsolve can be used effectively for this PDE system,

How to define the BCs correctly for the H-shaped cavity (especially the partially heated wall, θ=1 for 0.4 < y < 0.6), and

How to extract and plot the temperature (θ) and streamfunction (ψ) fields similar to the paper’s plots.

error_in_cavity_work.mw

How to solve this equation to get the the similar plots given in pdf peristatic_flow_work.pdf

here is the work sheet for the problem which is given in the pdf 

Download pulatile_flow_error.mw

I want the contour plots of the above worksheet

Good day.

I am looking into the behaviour of a function, V, that depends on several parameter values; these values are fixed in the attached example. However, I have encountered an issue that is puzzling me and I was hoping that someone may be able to shine a light on this for me.

Basically,  I would like to understand how the solution, V, behaves as the value of exponent, beta, approaches infinity.

Straightforward analysis suggests that the value of V tends to -10 as beta grows infinitely large (s -> C, and beta ->infinity and so, V -> -10 ... so far, so good).

However, when the function is plotted, the solution seems to converge to a value, 6.5, as beta tends to a very large number (10^17).

Now .. here's the mystery .. there appears to be a critical value of around 7.854 x 10^17; here, the limit seems to switch from V=6.5 to -10. Does this phenomenon correspond to a discontinuity or is it related to the computational process? Are there any built-in routines in Maple to check for such potential conditions?

Thanks for reading!

MaplePrimes_Sep_19.mw

dear sir here i am giving python code exicuting 3d plots but i cant plot animations like

restart;
Nc := .3; M := 1.5; QG := 1.5; Thetaa := .2; n1 := 1; n2 := 0; lambda1 := .1; lambda2 := .1; lambda3 := .1;

p := 2; a := .5; alpha1 := (1/2)*Pi;

p1 := 0.1e-1; p2 := 0.1e-1;
rf := 997.1; kf := .613; cpf := 4179; betaf := 21*10^(-5);

betas1 := .85*10^(-5); rs1 := 3970; ks1 := 40; cps1 := 765;
betas2 := 1.67*10^(-5); rs2 := 8933; ks2 := 401; cps2 := 385;

z1 := 1/((1-p1)^2.5*(1-p2)^2.5);
knf := kf*(ks1+2*kf-2*p1*(kf-ks1))/(ks1+2*kf+p1*(kf-ks1)); khnf := knf*(ks2+2*knf-2*p2*(knf-ks2))/(ks2+2*knf+p2*(knf-ks2));
z2 := 1-p1-p2+p1*rs1/rf+p2*rs2/rf;
z3 := 1-p1-p2+p1*rs1*cps1/(rf*cpf)+p2*rs2*cps2/(rf*cpf);
z4 := khnf/kf;
z5 := 1-p1-p2+p1*rs1*betas1/(rf*betaf)+p2*rs2*betas2/(rf*betaf);
OdeSys := z4*(X*(diff(Theta(X, tau), X, X))+diff(Theta(X, tau), X))/z3-(diff(Theta(X, tau), tau))-Nc*(Theta(X, tau)-Thetaa)^2*z5/z1-M^2*(Theta(X, tau)-Thetaa)^(p+1)*(1+n1*z2/z3)/(z3*(1-Thetaa)^p)-Nr*(Theta(X, tau)^4-Thetaa^4)/z3+QG*X*(1+lambda1*(Theta(X, tau)-Thetaa)+lambda2*(Theta(X, tau)-Thetaa)^2+lambda3*(Theta(X, tau)-Thetaa)^3)/z3; Cond := {Theta(0, tau) = 1+a*sin(alpha1), Theta(X, 0) = 0, (D[1](Theta))(1, tau) = 0};
OdeSys1 := z4*(X*(diff(Theta(X, tau), X, X))+diff(Theta(X, tau), X))/z3-(diff(Theta(X, tau), tau))-Nc*(Theta(X, tau)-Thetaa)^2*z5/z1-M^2*(Theta(X, tau)-Thetaa)^(p+1)*(1+n2*z2/z3)/(z3*(1-Thetaa)^p)-Nr*(Theta(X, tau)^4-Thetaa^4)/z3+QG*X*(1+lambda1*(Theta(X, tau)-Thetaa)+lambda2*(Theta(X, tau)-Thetaa)^2+lambda3*(Theta(X, tau)-Thetaa)^3)/z3; Cond1 := {Theta(0, tau) = 1+a*sin(alpha1), Theta(X, 0) = 0, (D[1](Theta))(1, tau) = 0};
OdeSysa := z4*(X*(diff(Theta(X, tau), X, X))+diff(Theta(X, tau), X))/z3-(diff(Theta(X, tau), tau))-Nc*(Theta(X, tau)-Thetaa)^2*z5/z1-M^2*(Theta(X, tau)-Thetaa)^(p+1)*(1+n1*z2/z3)/(z3*(1-Thetaa)^p)-Nr*(Theta(X, tau)^4-Thetaa^4)/z3+QG*X*(1+lambda1*(Theta(X, tau)-Thetaa)+lambda2*(Theta(X, tau)-Thetaa)^2+lambda3*(Theta(X, tau)-Thetaa)^3)/z3; Conda := {Theta(0, tau) = 1+a*cos(alpha1), Theta(X, 0) = 0, (D[1](Theta))(1, tau) = 0};
OdeSysa1 := z4*(X*(diff(Theta(X, tau), X, X))+diff(Theta(X, tau), X))/z3-(diff(Theta(X, tau), tau))-Nc*(Theta(X, tau)-Thetaa)^2*z5/z1-M^2*(Theta(X, tau)-Thetaa)^(p+1)*(1+n2*z2/z3)/(z3*(1-Thetaa)^p)-Nr*(Theta(X, tau)^4-Thetaa^4)/z3+QG*X*(1+lambda1*(Theta(X, tau)-Thetaa)+lambda2*(Theta(X, tau)-Thetaa)^2+lambda3*(Theta(X, tau)-Thetaa)^3)/z3; Conda1 := {Theta(0, tau) = 1+a*cos(alpha1), Theta(X, 0) = 0, (D[1](Theta))(1, tau) = 0};

colour := [cyan, black];
colour1 := [red, blue];
NrVals := [2.5, 3.5];
for j to numelems(NrVals) do Ans := pdsolve((eval([OdeSys, Cond], Nr = NrVals[j]))[], numeric, spacestep = 1/200, timestep = 1/100); Plots[j] := Ans:-plot(Theta(X, tau), tau = .8, linestyle = "solid", labels = ["Y", Theta(Y, tau)], color = colour[j], 'axes' = 'boxed'); eta[j] := Ans:-plot((int(z3*z5*Nc*(Theta(X, tau)-Thetaa)^2/(z1*z4)+NrVals[j]*(Theta(X, tau)^4-Thetaa^4)/z4+M^2*(Theta(X, tau)-Thetaa)^(p+1)*(1+n1*z2/z3)/(z4*(1-Thetaa)^p), X = 0 .. 1))/(z3*z5*Nc*(1-Thetaa)^2/(z1*z4)+NrVals[j]*(-Thetaa^4+1)/z4+M^2*(1-Thetaa)*(1+n1*z2/z3)/z4), tau = 0 .. 2, X = .8, linestyle = "solid", 'axes' = 'boxed', labels = [" &tau; ", " &eta; "], color = colour[j]); Ans1 := pdsolve((eval([OdeSys1, Cond1], Nr = NrVals[j]))[], numeric, spacestep = 1/200, timestep = 1/100); Plots1[j] := Ans1:-plot(Theta(X, tau), tau = .8, linestyle = "dash", labels = ["Y", Theta(Y, tau)], color = colour[j], 'axes' = 'boxed'); eta1[j] := Ans1:-plot((int(z3*z5*Nc*(Theta(X, tau)-Thetaa)^2/(z1*z4)+NrVals[j]*(Theta(X, tau)^4-Thetaa^4)/z4+M^2*(Theta(X, tau)-Thetaa)^(p+1)*(1+n2*z2/z3)/(z4*(1-Thetaa)^p), X = 0 .. 1))/(z3*z5*Nc*(1-Thetaa)^2/(z1*z4)+NrVals[j]*(-Thetaa^4+1)/z4+M^2*(1-Thetaa)*(1+n2*z2/z3)/z4), tau = 0 .. 2, X = .8, linestyle = "dash", 'axes' = 'boxed', labels = [" &tau; ", " &eta; "], color = colour[j]); Ansa := pdsolve((eval([OdeSysa, Conda], Nr = NrVals[j]))[], numeric, spacestep = 1/200, timestep = 1/100); Plotsa[j] := Ansa:-plot(Theta(X, tau), tau = .8, linestyle = "solid", labels = ["Y", Theta(Y, tau)], color = colour1[j], 'axes' = 'boxed'); etaa[j] := Ansa:-plot((int(z3*z5*Nc*(Theta(X, tau)-Thetaa)^2/(z1*z4)+NrVals[j]*(Theta(X, tau)^4-Thetaa^4)/z4+M^2*(Theta(X, tau)-Thetaa)^(p+1)*(1+n1*z2/z3)/(z4*(1-Thetaa)^p), X = 0 .. 1))/(z3*z5*Nc*(1-Thetaa)^2/(z1*z4)+NrVals[j]*(-Thetaa^4+1)/z4+M^2*(1-Thetaa)*(1+n1*z2/z3)/z4), tau = 0 .. 2, X = .8, linestyle = "solid", 'axes' = 'boxed', labels = [" &tau; ", " &eta; "], color = colour1[j]); Ansa1 := pdsolve((eval([OdeSysa1, Conda1], Nr = NrVals[j]))[], numeric, spacestep = 1/200, timestep = 1/100); Plotsa1[j] := Ansa1:-plot(Theta(X, tau), tau = .8, linestyle = "dash", labels = ["Y", Theta(Y, tau)], color = colour1[j], 'axes' = 'boxed'); etaa1[j] := Ansa1:-plot((int(z3*z5*Nc*(Theta(X, tau)-Thetaa)^2/(z1*z4)+NrVals[j]*(Theta(X, tau)^4-Thetaa^4)/z4+M^2*(Theta(X, tau)-Thetaa)^(p+1)*(1+n2*z2/z3)/(z4*(1-Thetaa)^p), X = 0 .. 1))/(z3*z5*Nc*(1-Thetaa)^2/(z1*z4)+NrVals[j]*(-Thetaa^4+1)/z4+M^2*(1-Thetaa)*(1+n2*z2/z3)/z4), tau = 0 .. 2, X = .8, linestyle = "dash", 'axes' = 'boxed', labels = [" &tau; ", " &eta; "], color = colour1[j]) end do;
plotA := plots:-display([seq(Plots[j], j = 1 .. 2)]);
plotB := plots:-display([seq(Plots1[j], j = 1 .. 2)]);
plotAA := plots:-display([seq(Plotsa[j], j = 1 .. 2)]);
plotBB := plots:-display([seq(Plotsa1[j], j = 1 .. 2)]);
plotC := plots:-display([seq(eta[j], j = 1 .. 2)]);
plotD := plots:-display([seq(eta1[j], j = 1 .. 2)]);
plotCC := plots:-display([seq(etaa[j], j = 1 .. 2)]);

plotDD := plots:-display([seq(etaa1[j], j = 1 .. 2)]);
plots:-display([plotA, plotB, plotAA, plotBB], size = [700, 700]);

plots:-display([plotC, plotD, plotCC, plotDD], size = [700, 700]);

how to take animation of the plots

i have seen some plots in maple also for that reason i have posted this question here

paper2_new_efficiency_plots_2025.mw